Making sessile drops easier 
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Abstract 



Using an identity, directly derived from tiie Young-Laplace equation, 
the problem of the equilibrium shape of an axisymmetric sessile drop 

rrt is reduced to a one-parameter shooting method problem. Based on the 

r^ method the numerical solutions for drops with Bond number up to 15 are 

*7^ plotted. The agreement between the method and the ADSA-D method 

'"rt as well as the experimental data is tested. A Mathematica code based on 

Pj the method is presented. 
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1 Introduction 

The problem of a drop on a horizontal surface with the effect of surface tension 
being balanced with gravity has been studied for more than a century. The 
early numerical solutions go back to 1883 yy, with updates by different authors 
[21 [3111] . Different perturbative treatments of the problem have been developed 
over the years, among them are those by [51 [HI [3 [HI H] for small drops (small Bond 
number). As large drops (or vanishing surface tension drops) are theoretically 
an infinitely large and thin film of liquid subjected to the boundary conditions at 
the outer edge, the limit of large Bond number falls and has been studied in the 
context of singular perturbation problems ^lOj . Based on the similarity between 
the truncated oblate spheroid and drop's shape, an approximated profile is 
suggested in [IT for the shape of the drop. In [T^] a new numerical treatment of 
the problem is given based on a variational method to minimize the total energy 
of the drop, by which the use of the tables by [1] is more direct than the earlier 
treatments. As another effort in this direction, in [13| the singular perturbation 
technique is used to obtain the asymptotic expressions describing the shape of 
small sessile and pendant drops. The study of the profiles of resting drops in 
different situations is particularly important for practical purposes. In fact, one 
of the most common methods to measure the surface tension of liquids is based 
on the matching between calculated drop's profiles and measured drop's shapes. 
Over the years, the optimization of matching methods between the calculated 
profiles and the experimental data on drop's profiles has been the subject of 
several research pieces [TJl [151 [T^ . 

Treating a sessile drop as a boundary value problem, it is shown that the 
multi-parameter shooting method is applicable (see e.g. [171 [181 [19] ) . It is the 
purpose of this note to show that a one-parameter shooting method is applicable 
to sessile drops as well. The reduction of parameters is the result of using an 
identity directly derived from the Young-Laplace equation |20| . By this identity 
the numerical procedure to reach the solution can be controlled by only one 
shooting parameter. 

The scheme of the rest of this note is as follows. In Sec. 2 the basic notions 
are shortly reviewed. In Sec. 3 the mathematical setup as well as the derivation 
of the mentioned identity is presented. Sec. 4 is devoted to some results and 
comparisons with some available numerical and experimental data. In Appendix 
a code in Mathematica based on the developed shooting procedure is presented. 

2 Basic notions 

The shape of a drop of liquid on a solid surface, in the idealized case (absence of 
impurities and pinning effects), is determined by the quantities: 1) the surface 
tension of liquid 7, 2) the solid-liquid adhesion coefficient a, 3) the shape of the 
solid surface, and due to the weight, 4) the drop's volume. At the solid-liquid- 
vapor (s.l.v.) point of contact, the contact angle i) in the equilibrium condition 
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Figure 1: Three classes of the drop's shape on a soUd surface. 



is given by the Young equation 
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cos?9=--l. (1) 
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Three classes of possibilities for the contact angle are presented in Fig. 1. 
At every point on the drop surface the Young-Laplace relation holds 

in which Ap = pi — p^ is the pressure jump across the surface, and (i?i,i?2) 
are two principal radii of curvature of the surface at the point. Provided by the 
hydrostatic laws, Ap is expressed in terms of the surface equation. Hence, the 
Young-Laplace relation is the partial differential equation which, accompanied 
by appropriate boundary conditions, determines the shape of the drop's surface. 
Heuristically, the final shape of drop is the result of the balance between the 
surface effects and the bulk ones. While the surface tension tends to decrease 
the surface of the drop, the adhesion coefficient tends to increase the surface of 
the contact region, and the gravity tends to lower the center of mass of the drop. 
In many practical cases the surface tension and the adhesion coefficient, though 
with opposite effects, may be considered at the same order, meaning 7 and a 
are comparable. For a drop with volume V and density g, the so-called Bond 
number defined by the dimensionless combination V'^'^gg/j would determine 
whether weight has the dominant contribution or not. When the weight is 
ignorable, only the contribution from the surface effects exist. Minimizing the 
area for a fixed volume, the drop's surface is part of a sphere. 

3 The mathematical setup 

Using the cylindrical coordinate setup given in Fig. 2, the total curvature of a 
surface with azimuthal symmetry, represented by z = f{p) is given by ^20 

where /' — df /dp. The pressure jump in presence of gravity gets contribution 
from the weight of the drop's layers as well, leading to 

Ap{z) = Ap^ + gg{h - z) (4) 
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Figure 2: The geometry of the mathematical setup for: a) i? < 90°, b) ?? > 90°. 



in which h is the height of the drop's apex, and Ap-y is a constant representing 
the pressure jump due to the surface tension. So, the Young-Laplace relation 
reads 
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(5) 



in which /+ and /_ are denoting the upper and lower parts of the drop, respec- 
tively (see Fig. 2b), and k :— Ap^/(27). At apex {h — /+(0)) we have i?i = i?2j 
and so by (l2|), K is simply the curvature at apex of drop. We mention /:^ < 
and f'_ > 0. The boundary conditions for {} > 90° are: 



/;(o) = 

f_ipo) = -tand, 
f-iPo) = 0. 



(6) 
(7) 
(8) 



In case with t? < 90° and ^ are valid for /+. 

The main issue with equation ([5]) is that the parameters h and k are not 
known at the first place, and would be determined only after the complete 
solution is available. So at starting point the main equation is not fully known. 
Further, the contact radius pa (Fig. 2), as the limiting value for the variable p, 
is not known at first place. These all indicate that the Young-Laplace equation 
can not be treated as straightforward as an ordinary boundary value problem. 
As will be seen shortly, by integrating the Young-Laplace equation an identity 
is obtained which relates the three unknown parameters in a very helpful way. 
In what follows we mainly consider the case with t? > 90°. The generalization 
to case with i) < 90° is rather straightforward. Integrating the Young-Laplace 
equation for the upper and lower parts of drop leads to [20] 
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in which we have used 
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l + f- 



— tan^ 
Vl + tan^ ^ 
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for I? > 90°. Subtracting ^ and ^ gives [5^ 
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in which we have used the relation for the volume of drop, 

TT = P/+(p)dp- / pf-{p)dp. 

27r Jo Jp„ 



(11) 



(12) 



(13) 



It is easy to show that identity ( 12 ) is valid for the acute contact angle (i? < 90°) 



as well. It is reminded that in obtaining (12 1 no approximation is used, and so 
it is an exact relation. 



Now, using the identity (12), the combination k and h in right-hand side of 
([5| can be replaced in favor of po, and the Young-Laplace equation turns to: 
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In above the only unknown parameter is po. For the case with d < 90° only /+ 
in above should be kept. 

For later use, let us remind the spherical solution of weightless drop. It is 
easy to check that, by setting g = and k = 1/R, the expression [20] 



2 = fo±{p) = zo± VR^ - P^i 
satisfies ((sl) . The radius R is fixed by the volume of spherical cap, 

V ^ ^i?^(l-cosi?)2(2 + cos^). 

Following a simple geometrical argument in the sphere (Fig. 2), we have 



Po ~ Rsin'd 

zq = —Rcosd 

Pi = R 

h = R + Zq 



3 = 



(15) 
(16) 

(17) 



As in case with i) > 90° around the equatorial radius pi (Fig. 2b) /j_ —>■ =F0O: 
it would be convenient to switch to the spherical coordinates in which the whole 
surface of drop is covered by one function. Among many others, we found the 
choice illustrated in Fig. 3 more convenient, in which 



2 = r{9) sin 6*, 
p — r{6) cos 9. 



(18) 
(19) 
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Figure 3: The geometry with the polar coordinate (r, 
the sphere solution of weightless drop. 



The dashed curve is 



In the new coordinates the Young-Laplace equation takes the following form 



r r' cos 9 
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(20) 



in which r' = dr/d9. Using the identity ( 12 ) leads to 
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(21) 



in which c — gQ/"f, the so-called the capillary constant. As we are considering 
axisymmetric drops, it is sufficient to restrict the polar coordinate 9 only to the 
left-side of drop, 



0< 6l< 90°. 



(22) 



Using ([BJI-Q and (18)-(19) the boundary conditions in new coordinates are 
found 



r'(0) = -pocot??, 
riO) = Po, 
r'(90°) = 0. 



(23) 
(24) 
(25) 



Therefore the shooting method b y on e controlling parameter po can be employed 
to obtain numerical solutions of (21 1. The value of pp" = i?sin'i9 of weightless 
drop can be chosen as the initial guess of the shooting parameter po- As it is 
known that, due to weight, the radius of contact region increases (see e.g. pUj). 
upon stepwise increasing of the parameter po smd checking the condition ( |25[ ) 
at apex of drop, one can reach the desired accuracy. 




Figure 4: The plots of solutions of ( |2l[ ) for drops with acute contact angle but 
different Bond numbers (given on each plot). In all samples: V = 1, 7 = 70, 
g = I, and 'd — 45° (the varying input value is g). Scales on axes: 1:1. 




Figure 5: The same as Fig. 4 but with obtuse contact angle 1} = 145° 



4 Samples of solutions and comparisons 

In order to assess the accuracy of the present method, the results based on it 
are compared with some available data. However, it would be illustrative to 
begin with a demonstration of the results based on the method. The numerical 
solutions of (21 1 with different Bond numbers for 1? = 45° and d = 145° are 
presented in Fig. 4 and Fig. 5, respectively. The presented plots cover Bond 
numbers from to 15. As expected, by increasing the Bond number, the apex's 
height decreases while contact radius, as well as equatorial radius for case with 
1} > 90°, increase. 

As a further test, in Tab. 1 the results of application of one-parameter shoot- 
ing method to (21) are compared with those by p2] by the ADSA-D method. 



The input values used by T^ are taken from the data generated by ALFI pro- 
gram. It is mentioned that the present method returns the correct values with 
negligible errors. 

The comparison with some available experimental data with the values by 
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c 


V 
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flmax. 


K, 


flmax 




cm^^ 


em"^ 


em~^ 


cm 


(err.) 


(err.) 


Sl/120° 


13.45 


0.4578 


0.7340 


0.6562 


0.7338 
(0.03%) 


0.6562 
(0.00%) 


S2/10° 


13.45 


0.1236 


0.030 


1.1152 


0.030 
(0.00%) 


1.1150 

(0.02%) 


S3/80° 


13.45 


0.0014 


10.0 


0.0965 


10.09 

(0.90%) 


0.0957 

(0.84%) 


S4/30° 


1.00 


0.0486 


1.00 


0.4849 


1.00 
(0.00%) 


0.4846 
(0.06%) 


S5/45° 


1000 


0.00089 


1.00 


0.1299 


1.00 


0.1299 



(0.00%) (0.00%) 

Table 1: The comparison between the numerical solutions by [18] and those by the 
present work. The values by [18] are taken from Tabs. 1 & 2. The parameter c, as in 
this work, denotes the capillary constant gg/'^. The curvature at apex (k) is denoted 
in [18] by h. i?max. denotes the contact radius po for •& < 90° and the equatorial radius 
pi for 1? > 90° (Fig. 2). 



solutions of (21 1 are presented in Tab. 2 and Tab. 3 for acute and obtuse contact 



angles, respectively. 
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0.253 


0.454 
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0.245 
(3.2%) 



Table 2: The comparison between experimental data for drops with i? < 90° 



and values by numerical solutions of (211. For water drops: g = 0.997 g/cm , 
7 = 72.0 mj/m . For formamide drops: g — 1.133 g/cm , 7 = 58.2 mj/m . For 
all drops: g — 980.7 cm/s . 



A lyiatliematica code 

Here a code in Mathematica based on the method is presented. In the below the 
following input values from Tab. 3 are taken: V = 0.0892 cm'^, 7 — 72.0 mJ/m , 
g = 980.7 cm/s , g = 0.997 g/cm , d = 117.34°. The increasing step in po is 
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0.1157 


0.1435 


0.1950 


0.1160 
(0.3%) 


0.1414 
(1.5%) 


0.1930 

(1.0%) 



Table 3: The comparison between experimental data for drops with -3 > 90° and vab 
ues by numerical solutions of (21 1. For water drops: g = 0.997 g/cm , 7 = 72.0 mj/m . 
For mercury drops: g — 13.55 g/cm , 7 = 480 mJ/m . For all drops: g — 980.7 cm/s . 



taken equal to 0.00001, with the criteria that slope at apex would be less than 
0.1. The outputs of the code, together with the function r{6) (ryl[tha]), are: 
the resulting slope at the apex (slope), the contact radius po (^hO), the apex's 
height (h), the angel di, radius pi, and height hi at equator (thl, rhl and hi), 
together with the plot of the drop. 



V=0 . 0892 ; gaimna=72 ; g=980 . 7 ; dens=0 . 997 ; vth=117 . 34 Degree ; 

R=(3 V/(Pi(l-Cos[vth])^2(2+Cos[vth])))^(l/3) ; 

c=g dens/gamma ;rhOsph=R Sin[vth]; 

thinin=0 Degree ;thmax=89 .99 Degree; 

ylsol[rhO_] :=NDSolve [{D[r [th] r' [th] Cos [th] /Sqrt [r [th]^2+r ' [th]^2] ,th] 

-(2 r[th]'^2+r' [th] ^2)Cos [th] /Sqrt [r [th]^2+r ' [th]^2] 

+2(Sin[vth]/rhO+c V/(2 Pi rhO'^2) )r [th] ^^2 Cos[th] 

-c r[th]^3Sin[th] Cos[th]==0, 

r [thmin] ==rhO , r ' [thmin] ==-rhO Cot [vth] } , r , {th , thmin , thmax}] ; 

slope [rhO_] : =r ' [thmax] / . ylsol [rhO] [ [1] ] //Quiet ; 

rhOt=rhOsph ; While [Abs [slope [rhOt] ] >0.1 , rhOt=rhOt+0 . 00001] ; 

ryl [tha J : =r [tha] / . ylsol [rhOt] [ [1] ] ; 

Print ["slope=" , slope [rhOt] //N] 

apx=ryl [thmax] //N ; 

Print ["rhO=",rhOt," h=",apx] 

thl=tha/ . FindRoot [D [ryl [tha] Cos [tha] ,tha] , {tha, 0.05}] [[1]] ; 

rhl=ryl [thl] Cos [thl] //N ; 

hl=ryl [thl] Sin [thl] //N ; 

Print ["thl=", thl," rhl=",rhl," hl=",hl] 

ParametricPlot [{ryl [tha] Cos [tha] ,ryl [tha] Sin [tha] } , {tha, thmin, thmax}, 

AspectRatio— >Automatic,PlotStyle— )-Black] 
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